Nrf2 BL6 PEITC 16S Microbiome Data Visualization

Date: Mon Sep 2 11:47:39 2019
Scientist: Ran Yin
Sequencing (Waksman): Dibyendu Kumar
Statistics: Davit Sargsyan
Principal Investigator: Ah-Ng Kong

1 Meta data

options(stringsAsFactors = FALSE,
        scipen = 999)
# # Increase mmemory size to 64 Gb----
# invisible(utils::memory.limit(65536))
# str(knitr::opts_chunk$get())
# # NOTE: the below does not work!
# knitr::opts_chunk$set(echo = FALSE, 
#                       message = FALSE,
#                       warning = FALSE,
#                       error = FALSE)
require(knitr)
require(kableExtra)
require(data.table)
require(phyloseq)
require(ggplot2)
require(plotly)
require(DT)
require(shiny)
source("source/functions_may2019.R")
# On Windows set multithread=FALSE----
mt <- TRUE

# Load data----
# Counts
load("data_may2019/ps_may2019.RData")
# Taxonomy
load("data_may2019/taxa.RData")
taxa <- data.table(seq16s = rownames(taxa),
                   taxa)
# Samples
samples <- ps_may2019@sam_data
datatable(samples,
              options = list(pageLength = nrow(samples)))

# Taxonomic Ranks:
# **K**ing **P**hillip **C**an n**O**t **F**ind **G**reen **S**ocks
# * Kingdom                
# * Phylum                    
# * Class                   
# * Order                   
# * Family     
# * Genus     
# * Species  

2 Prune data

Keep only bacterial OTUs.

Check bacterial distribution by phylum.

First, all OTUs that were not bacteria were removed. Also, OTUs that were not mapped to a phylum were removed as these are usually sequencing artifacts.

3 Richness (Alpha diversity)

ps0@sam_data$Diet_Week <- paste(samples$TREATMENT,
                                samples$WEEK,
                                sep = "_")
ps0@sam_data$ID <- substr(x = ps0@sam_data$SAMPLE_NAME,
                          start = 2,
                          stop = 3)
p1 <- plot_richness(ps0,
                    x = "Diet_Week", 
                    measures = "Shannon") +
  geom_point(aes(fill = ID),
             shape = 21,
             size = 3,
             color = "black") +
  geom_line(aes(group = ID),
            color = "black") +
  scale_x_discrete("") +
  theme(axis.text.x = element_text(angle = 45))
ggplotly(p = p1,
         tooltip = c("ID",
                     "value"))

4 OTU table (first 10 rows)

5 Total counts per sample (i.e. sequencing depth)

6 Counts at Phylum level

7 Relative abundance (%) at Phylum level

Remove phyla with relative abundance of >= 1% in less than 10% of samples.

t1 <- data.table(Phylum = ra_p$Phylum,
                 `Number of Samples` = rowSums(ra_p[, 2:ncol(ra_p)] >= 0.01))
t1$`Percent Samples` <-  t1$`Number of Samples`/30
setorder(t1, -`Number of Samples`)
datatable(t1,
          rownames = FALSE,
          caption = "Taxonomic  count table",
          class = "cell-border stripe",
          options = list(search = FALSE,
                         pageLength = nrow(t1))) %>%
  formatPercentage(columns = 3,
                   digits = 1)

Hence, only 6 out of 13 Phyla were studied in this analysis: Actinobacteria, Bacteroidetes, Firmicutes, Proteobacteria, Tenericutes and Verrucomicrobia. Relative abundance at the next taxonomic level (Class) was, therefore, computed relative to the sum of these 6 Phyla.

[1] "Bacteroidetes, Firmicutes, Proteobacteria, Verrucomicrobia, Actinobacteria, Tenericutes"

7,628 OTUs, down from 7,764 OTUs in the previous table.

8 Relative Abundance in Samples at Different Taxonomic Ranks

8.1 1. Class

p0 <- ggplot(mu,
             aes(x = Week,
                 y = x,
                 group = Treatment)) +
  facet_wrap(~ rev(Class),
             scale = "free_y") +
  geom_line() +
  geom_point(aes(shape = Treatment,
                 color = Treatment),
             size = 5,
             alpha = 0.5) +
  scale_y_continuous("Relative abundance (%)")
print(p0)

p1 <- ggplot(mu,
             aes(x = x,
                 y = Class,
                 fill = Treatment,
                 shape = Week)) +
  # facet_wrap(~ Sex, nrow = 1) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)")
ggplotly(p1)

8.2 2. Order

p0 <- ggplot(mu,
             aes(x = Week,
                 y = x,
                 group = Treatment)) +
  facet_wrap(~ Order,
             scale = "free_y") +
  geom_line() +
  geom_point(aes(shape = Treatment,
                 color = Treatment),
             size = 5,
             alpha = 0.5)
print(p0)

p1 <- ggplot(mu,
             aes(x = x,
                 y = Order,
                 fill = Treatment,
                 shape = Week)) +
  # facet_wrap(~ Sex, nrow = 1) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)")
ggplotly(p1)

9 Session Information

sessionInfo()

R version 3.5.0 (2018-04-23)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.5 (Maipo)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] shiny_1.3.2       DT_0.6            plotly_4.9.0      ggplot2_3.2.0     phyloseq_1.26.1   data.table_1.12.2 kableExtra_1.1.0 
[8] knitr_1.23       

loaded via a namespace (and not attached):
 [1] Biobase_2.42.0      httr_1.4.0          tidyr_0.8.3         jsonlite_1.6        viridisLite_0.3.0   splines_3.5.0      
 [7] foreach_1.4.4       assertthat_0.2.1    stats4_3.5.0        yaml_2.2.0          pillar_1.4.2        lattice_0.20-35    
[13] glue_1.3.1          digest_0.6.19       promises_1.0.1      XVector_0.22.0      rvest_0.3.4         colorspace_1.4-1   
[19] httpuv_1.5.1        htmltools_0.3.6     Matrix_1.2-14       plyr_1.8.4          pkgconfig_2.0.2     zlibbioc_1.28.0    
[25] xtable_1.8-4        purrr_0.3.1         scales_1.0.0        webshot_0.5.1       later_0.8.0         tibble_2.1.3       
[31] mgcv_1.8-23         IRanges_2.16.0      withr_2.1.2         BiocGenerics_0.28.0 lazyeval_0.2.2      mime_0.6           
[37] survival_2.41-3     magrittr_1.5        crayon_1.3.4        evaluate_0.14       nlme_3.1-137        MASS_7.3-49        
[43] xml2_1.2.0          vegan_2.5-5         tools_3.5.0         hms_0.4.2           stringr_1.4.0       Rhdf5lib_1.4.3     
[49] S4Vectors_0.20.1    munsell_0.5.0       cluster_2.0.7-1     Biostrings_2.50.2   ade4_1.7-13         compiler_3.5.0     
[55] rlang_0.4.0         rhdf5_2.26.2        grid_3.5.0          iterators_1.0.10    biomformat_1.10.1   rstudioapi_0.10    
[61] htmlwidgets_1.3     igraph_1.2.4.1      rmarkdown_1.13      gtable_0.3.0        codetools_0.2-15    multtest_2.38.0    
[67] reshape2_1.4.3      R6_2.4.0            dplyr_0.8.0.1       permute_0.9-5       readr_1.3.1         ape_5.3            
[73] stringi_1.4.3       parallel_3.5.0      Rcpp_1.0.0          tidyselect_0.2.5    xfun_0.7